// #include "FuncFiller.H"
// #include <array>
#include "Circle.H"
// #include "EmbededGrid.H"
// #include "DirectPossionSolver.H"

enum {Dim = 2};
using rVec = Vec<Real, Dim>;
using iVec = Vec<int, Dim>;


// Real PHI(const rVec &x)
// {
//    return sin(4*M_PI*x[0]) * cos(3*M_PI*x[1]);
// }

// rVec DPHI(const rVec &x)
// {
//   return rVec{
//       4 * M_PI * cos(4 * M_PI * x[0]) * cos(3 * M_PI * x[1]),
//       -3 * M_PI * sin(4 * M_PI * x[0]) * sin(3 * M_PI * x[1])};
// }

// Real DDPHI(const rVec &x)
// {
//   return -25 * M_PI * M_PI * sin(4 * M_PI * x[0]) * cos(3 * M_PI * x[1]);
// }

// Real negDDPHI(const rVec &x)
// {
//   return 25 * M_PI * M_PI * sin(4 * M_PI * x[0]) * cos(3 * M_PI * x[1]);
// }

// Real BETA(const rVec &x)
// {
//   return 1.0 + x[0] + x[1];
// }

// rVec DBETA(const rVec &x)
// {
//   return rVec{1.0, 1.0};
// }

// Real DDBETA(const rVec &x)
// {
//   return 0;
// }

// Real EXPT1(const rVec &x)
// {
//   return exp(sin(x[0])+x[1]);
// }

// rVec DEXPT1(const rVec &x)
// {
//   return rVec{cos(x[0])*exp(sin(x[0])+x[1]), exp(sin(x[0])+x[1])};
// }

// Real negDDEXPT1(const rVec &x)
// {
//   return -exp(sin(x[0])+x[1])*(1.0+cos(x[0])*cos(x[0])-sin(x[0]));
// }





void test1()
{
  using namespace std;
  OrientedCircle C1({0.3, 0.2}, 0.1,Clockwise);
  cout << C1.inEnclosure({0.3, 0.2}) << '\n';
  cout << C1.inEnclosure({0.3, 1e-5}) << '\n';
  cout << C1.inEnclosure({-5, -5}) << '\n';
  return;
//   vector<rVec> square1 = C1.minimalEmbededdSquare();
//   cout << square1[0] << square1[1] << '\n';
//   // cout<<C1.polarRep(M_PI)<<'\n';
//   cout << C1.getNormalGrad(C1.polarRep(M_PI / 2)) << '\n';
//   vector<rVec> cutPoint = C1.intersection2CartesianLine(0, 0.3);
//   cout << cutPoint[0] << cutPoint[1] << '\n';
//   cout << C1.intersection2CartesianLine(1, {0.21, 0.2}) << '\n';
//   Real rad1 = C1.getRadian({0.2,0.2});
//   cout<<rad1<<C1.polarRep(rad1);
};
// void test2()
// { /// Tensor支持负索引且不是越界的那种
//   Box<2> bx(-1, 2);
//   Tensor<Real, Dim> m(bx);
//   loop_box_2(bx, i, j)
//       m(i, j) = i + j;
//   std::cout << m << '\n';
// };
// void test3()
// {
//   /// Regular
//   int N = 8;
//   Box<Dim> b1(0, N);
//   RectDomain<Dim> pGrid(Box<Dim>(0, N), 1.0 / N);
//   ExcludedCircle<2> rce1(pGrid, (Circle){{0.5, 0.5}, 0.2}, 1e-8);

//   FuncFiller<Dim> ff(pGrid);

//   Box<Dim> gbx = pGrid;
//   Tensor<Real, Dim> phi(gbx);
//   Tensor<Real, Dim> err(gbx);
//   Tensor<Real, Dim> res(gbx);
//   Tensor<Real, Dim - 1> bcData[2 * Dim];

//   loop_box_2(gbx,i,j)
//   {
//     ff.fillDomainVal(phi,DDBETA);
//   }

//   for (int d = 0; d < Dim; ++d)
//   {
//     for (int side : {-1, 1})
//     {
//       int a = (side == -1) ? (gbx.lo()[d]) : (gbx.hi()[d]);
//       int m = 2 * d + (side + 1) / 2;
//       bcData[m].resize(reduce(gbx, m / 2));
//       ff.fillRectBdrVal(bcData[m], d, side, BETA);
//     }
//   }
//   DirectPossionSolver<2> DPS1(rce1, "DDDDD",bcData,phi,BETA);
// };
// void test4()
// {
//   int N = 8;
//   Box<Dim> b1(0, N);
//   RectDomain<Dim> pGrid(Box<Dim>(0, N), 1.0 / N);
//   FuncFiller<Dim> ff(pGrid);

//   Box<Dim> gbx = pGrid;
//   Tensor<Real, Dim> phi(gbx);
//   Tensor<Real, Dim> err(gbx);
//   Tensor<Real, Dim> res(gbx);
//   Tensor<Real, Dim - 1> bcData[2 * Dim];

//   for (int d = 0; d < Dim; ++d)
//   {
//     for (int side : {-1, 1})
//     {
//       int a = (side == -1) ? (gbx.lo()[d]) : (gbx.hi()[d]);
//       int m = 2 * d + (side + 1) / 2;
//       bcData[m].resize(reduce(gbx, m / 2));
//       ff.fillRectBdrVal(bcData[m], d, side, BETA);
//     }
//   }
//   for(int i:{0,1,2,3})
//   std::cout<<bcData[i]<<'\n';
// }
// void test5()
// {
//   /// IrRegular Simple Case
//   int N = 16;
//   Box<Dim> b1(0, N);
//   RectDomain<Dim> pGrid(Box<Dim>(0, N), 1.0 / N);
//   ExcludedCircle<2> rce1(pGrid, (OrientedCircle){{0.4, 0.5}, 0.2}, 1e-8);
//   FuncFiller<Dim> ff(pGrid);

//   Box<Dim> gbx = pGrid;
//   Tensor<Real, Dim> phi(gbx);
//   Tensor<Real, Dim> err(gbx);
//   Tensor<Real, Dim> res(gbx);
//   Tensor<Real, Dim - 1> bcData[2 * Dim];

//   loop_box_2(gbx,i,j)
//   {
//     ff.fillDomainVal(phi,DDBETA);
//   }

//   for (int d = 0; d < Dim; ++d)
//   {
//     for (int side : {-1, 1})
//     {
//       int a = (side == -1) ? (gbx.lo()[d]) : (gbx.hi()[d]);
//       int m = 2 * d + (side + 1) / 2;
//       bcData[m].resize(reduce(gbx, m / 2));
//       ff.fillRectBdrVal(bcData[m], d, side, BETA);
//     }
//   }
//   //ff.fillRectBdrNormalGrad(bcData[3],1,1,DBETA);
//   // ff.fillRectBdrVal(bcData[2], 1, -1, NGBETA2);
//   // ff.fillRectBdrVal(bcData[1], 0, 1, NGBETA3);
//  DirectPossionSolver<2> DPS1(rce1, "DDDDN",bcData,phi,[&](const auto &x){return dot(DBETA(x),rce1.embeddingCircle.getNormalGrad(x));});
//   std::cout<<DPS1.computeErrNorm(BETA,0)<<'\n';
//  DPS1.plotIn("file.m",BETA);
// };

// void test6()
// {
//   /// IrRegular
//   int N = 8;
//   std::cin>>N;
//   Box<Dim> b1(0, N);
//   RectDomain<Dim> pGrid(Box<Dim>(0, N), 1.0 / N);
//   ExcludedCircle<2> rce1(pGrid, (Circle){{0.4, 0.5}, 0.2}, 1e-8);
//   FuncFiller<Dim> ff(pGrid);
//   Box<Dim> gbx = pGrid;
//   Tensor<Real, Dim> phi(gbx);
//   Tensor<Real, Dim> err(gbx);
//   Tensor<Real, Dim> res(gbx);
//   Tensor<Real, Dim - 1> bcData[2 * Dim];
//   loop_box_2(gbx,i,j)
//   {
//     ff.fillDomainVal(phi,negDDPHI);
//   }
//   for (int d = 0; d < Dim; ++d)
//   {
//     for (int side : {-1, 1})
//     {
//       int a = (side == -1) ? (gbx.lo()[d]) : (gbx.hi()[d]);
//       int m = 2 * d + (side + 1) / 2;
//       bcData[m].resize(reduce(gbx, m / 2));
//       ff.fillRectBdrVal(bcData[m], d, side, PHI);
//     }
//   }

//   // ff.fillRectBdrNormalGrad(bcData[3],1,1,DPHI);
//   // ff.fillRectBdrNormalGrad(bcData[2], 1, -1, DPHI); 
//  //DirectPossionSolver<2> DPS1(rce1, "DDNND",bcData,phi,PHI);
//   DirectPossionSolver<2> DPS1(rce1, "DDDDN",bcData,phi,[&](const auto &x){return dot(DPHI(x),rce1.embeddingCircle.getNormalGrad(x));});
//  std::cout<<DPS1.computeErrNorm(PHI,0)<<'\n';
//  DPS1.plotIn("file6.m",PHI);
// };

// void test7()
// {/// tan(pi)=?
// Real rad = 0;
// rVec dx = {0.5,0.5};
//       Real l = dx[0]/std::max(abs(sin(rad)),abs(cos(rad)));
//       Real cutRadio = std::min(abs(tan(rad)),1.0/abs(tan(rad))); 
// std::cout<<l<<" "<<cutRadio<<'\n';
// }

// void test9()
// {
//   /// IrRegular
//   int N = 8;
//   std::cin>>N;
//   Box<Dim> b1(0, N);
//   RectDomain<Dim> pGrid(Box<Dim>(0, N), 1.0 / N);
//   ExcludedCircle<2> rce1(pGrid, (OrientedCircle){{0.4, 0.5}, 0.2}, 1e-8);
//   FuncFiller<Dim> ff(pGrid);
//   Box<Dim> gbx = pGrid;
//   Tensor<Real, Dim> phi(gbx);
//   Tensor<Real, Dim> err(gbx);
//   Tensor<Real, Dim> res(gbx);
//   Tensor<Real, Dim - 1> bcData[2 * Dim];
//   loop_box_2(gbx,i,j)
//   {
//     ff.fillDomainVal(phi,negDDEXPT1);
//   }
//   for (int d = 0; d < Dim; ++d)
//   {
//     for (int side : {-1, 1})
//     {
//       int a = (side == -1) ? (gbx.lo()[d]) : (gbx.hi()[d]);
//       int m = 2 * d + (side + 1) / 2;
//       bcData[m].resize(reduce(gbx, m / 2));
//       ff.fillRectBdrVal(bcData[m], d, side, EXPT1);
//     }
//   }

//   // ff.fillRectBdrNormalGrad(bcData[3],1,1,DEXPT1);
//   // ff.fillRectBdrNormalGrad(bcData[2], 1, -1, DEXPT1); 
//   // DirectPossionSolver<2> DPS1(rce1, "DDNND",bcData,phi,EXPT1);
//  DirectPossionSolver<2> DPS1(rce1, "DDDDN",bcData,phi,[&](const auto &x){return dot(DEXPT1(x),rce1.embeddingCircle.getNormalGrad(x));});
//  std::cout<<DPS1.computeErrNorm(EXPT1,0)<<'\n';
//  DPS1.plotIn("file9.m",EXPT1);
// };

// void test8()
// {
//   /// IrRegular
//   int N = 8;
//   std::cin>>N;
//   Box<Dim> b1(0, N);
//   RectDomain<Dim> pGrid(Box<Dim>(0, N), 1.0 / N);
//   ExcludedCircle<2> rce1(pGrid, (OrientedCircle){{0.4, 0.5}, 0.2}, 1e-8);
//   FuncFiller<Dim> ff(pGrid);
//   Box<Dim> gbx = pGrid;
//   Tensor<Real, Dim> phi(gbx);
//   Tensor<Real, Dim> err(gbx);
//   Tensor<Real, Dim> res(gbx);
//   Tensor<Real, Dim - 1> bcData[2 * Dim];
//   loop_box_2(gbx,i,j)
//   {
//     ff.fillDomainVal(phi,negDDPHI);
//   }
//   for (int d = 0; d < Dim; ++d)
//   {
//     for (int side : {-1, 1})
//     {
//       int a = (side == -1) ? (gbx.lo()[d]) : (gbx.hi()[d]);
//       int m = 2 * d + (side + 1) / 2;
//       bcData[m].resize(reduce(gbx, m / 2));
//       ff.fillRectBdrVal(bcData[m], d, side, PHI);
//     }
//   }
//  DirectPossionSolver<2> DPS1(rce1, "DDDDN",bcData,phi,[&](const auto &x){return dot(DPHI(x),rce1.embeddingCircle.getNormalGrad(x));});
//  std::cout<<DPS1.computeErrNorm(PHI,0)<<'\n';
//  DPS1.plotIn("file8.m",PHI);
// };

// template<class negDDExpr,class Expr,class DExpr>
// void testExample(const negDDExpr &negDDexpr,const Expr &expr,const DExpr Dexpr,const char *aBcTypes, int N = 8)
// {
//   char bcTypes[10];
//   strncpy(bcTypes, aBcTypes, 10);
//   Box<Dim> b1(0, N);
//   RectDomain<Dim> pGrid(Box<Dim>(0, N), 1.0 / N);
//   ExcludedCircle<2> rce1(pGrid, (OrientedCircle){{0.4, 0.5}, 0.2}, 1e-8);
//   FuncFiller<Dim> ff(pGrid);
//   Box<Dim> gbx = pGrid;
//   Tensor<Real, Dim> phi(gbx);
//   Tensor<Real, Dim - 1> bcData[2 * Dim];
//   loop_box_2(gbx,i,j)
//   {
//     ff.fillDomainVal(phi,negDDexpr);
//   }
//   for (int d = 0; d < Dim; ++d)
//   {
//     for (int side : {-1, 1})
//     {
//       int m = 2 * d + (side + 1) / 2;
//       switch (bcTypes[m])
//       {
//       case 'D':
//       {
//           int a = (side == -1) ? (gbx.lo()[d]) : (gbx.hi()[d]);
//           bcData[m].resize(reduce(gbx, m / 2));
//           ff.fillRectBdrVal(bcData[m], d, side, expr);
//         break;
//       }
//       case 'N':
//       {
//           ff.fillRectBdrNormalGrad(bcData[m],d,side,Dexpr);
//       }
//         break;
//       default:
//         break;
//       }
//     }
//   }
//   switch (bcTypes[4])
//   {
//   case 'D':
//   {
//     DirectPossionSolver<2> DPS1(rce1, bcTypes,bcData,phi,expr);  
//     std::cout<<DPS1.computeErrNorm(expr,0)<<'\n';
//     DPS1.plotIn("fileExample.m",expr);      
//     break;
//   }
//   case 'N':{
//     DirectPossionSolver<2> DPS1(rce1, bcTypes,bcData,phi,[&](const auto &x){return dot(Dexpr(x),rce1.embeddingCircle.getNormalGrad(x));});
//     std::cout<<DPS1.computeErrNorm(expr,0)<<'\n';
//     DPS1.plotIn("fileExample.m",expr);
//     break;
//   }
//   default:
//     break;
//   }
// }

int main()
{
    test1();
//   for(int i:{8,16,32,64})
//   testExample(negDDEXPT1,EXPT1,DEXPT1,"DDDDN",i);
}



